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Abstract 

We discuss how we formulate time evolution of physical quantities in the framework of the 
Rigged QED (Quantum Electrodynamics). The Rigged QED is a theory which has been proposed 
to treat dynamics of electrons, photons and atomic nuclei in atomic and molecular systems in 
a quantum field theoretic way. To solve the dynamics in the Rigged QED, we need different 
techniques from those developed for the conventional QED. As a first step toward this issue, 
we propose a procedure to expand the Dirac field operator, which represents electrons, by the 
electron annihilation/creation operators and solutions of the Dirac equation for electrons in nuclear 
potential. Similarly, the Schrodinger field operators, which represent atomic nuclei, are expanded 
by nucleus annihilation/creation operators. Then we derive time evolution equations for these 
annihilation and creation operators and discuss how time evolution of the operators for physical 
quantities can be calculated. In the end, we propose a method to approximate the evolution 
equations of the operators by the evolution equations for the density matrices of electrons and 
atomic nuclei. Under this approximation, we carry out numerical simulation of the time evolution 
of electron charge density of a hydrogen atom. 
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I. INTRODUCTION 



Recently, the technology for observation and manipulation of quantum systems involving 
the interaction between light and matter is developed with increasing speed and new ex- 
perimental tools are becoming available. Remarkable features of such developments include 

fin nn 

real-time observation [l|, |2j, single-photon generation [3|-[6[ and so on. Then, as for the theo- 
retical side, it would be desirable to develop a formalism to simulate dynamical phenomena 
observed by such experiments treating light and matter and the interaction between them 
as accurately as possible. That is, by Quantum Electrodynamics (QED). 

QED is probably the most successful fundamental physical theory we have. It is the 
theory of electrically charged particles and photons and interactions between them, taking 
the form of quantum theory of fields. Its accurate predictions include the Lamb shift of the 
energy levels of the hydrogen atom and the anomalous magnetic moment of the electron. 
In addition to such stationary properties, it can also predict dynamical properties like the 
cross section for the electron-positron scattering (Bhabha scattering) to a great accuracy. 
The successes are due to the development of the method to compute the QED interaction 
as a perturbation. The success of the perturbative approach, in turn, owes much to the 
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capability of preparing the asymptotic states as the unperturbed states (and o: 
the renormalizability) . It has been long since such formalism was established [7, 

However, use of the asymptotic states and subsequently the S"-matrix formulation has 
shortcomings when we would like to know the step-by-step time evolution of the quantum 
system employing QED. This is because it just describes the transition between infinite 
past ("in-state") and infinite future ("out-state"). Although it works fine for calculations of 
cross sections of scattering processes or energy shifts of bound states, it cannot be used for 
simulating the time evolution of the quantum system. 

Hence, so far, simulations of time evolution of the quantum system involving light and 
matter are performed using the time-dependent Schrod inge r or Dirac equation with classical 
electromagnetic fields (i.e. semi-classical treatment) 9Mll| or using a quantized photon field 



with the very much simplified matter part (so that the interaction between the photon and 
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14|. It is true that these approx- 



matter is introduced somewhat in an ad hoc manner) 
imations are appropriate for a wide range of systems. For example, the former treatment 
is reasonable for a photon field produced by laser and the latter treatment for a system in 
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cavity QED. We consider, however, it is of great importance to develop a simulation method 
based on QED in the form of quantum field theory, and without recourse to the perturbative 
approach. 1 

Besides the non-perturbative time evolution scheme, what we wish to add to the standard 
QED framework is the atomic nucleus. This is because, contrary to the high-energy physics 
application, the existence of the atomic nuclei is essential for atomic and molecular science 
in which we are interested. We would like to have a field theoretic QED formalism which is 
applicable to low energy phenomena such as chemical reactions and photoionization process. 



In this paper, as has been proposed in Ref. 
of freedom in the framework of the Rigged QED. 



181 ]. we include the atomic nucleus degree 



The Rigged QED 17H19| is a theory which has been proposed to treat dynamics of 
charged particles and photons in atomic and molecular systems in a quantum field theoretic 
way. In addition to the ordinary QED which is Lorentz invariant with Dirac (electron) 
field and U(l) gauge (photon) field, Schrodinger fields which represent atomic nuclei are 
added. In this way, dynamics of nuclei and their interaction with photons can be treated 
in a unified manner. Incidentally, we here note on semantics of the word "rigged" . It has a 
nautical connotation as in the phrase "fully rigged ship" and means "equipped" ? Namely, 
the Rigged QED is "QED equipped with Schrodinger fields" as we just described above. 

Including the atomic nucleus degrees of freedom as quantum fields is crucial. Since we 
would like to treat the interaction between an electron and an atomic nuclei or between 
two nuclei as the one which is mediated by the quantized photon field (not by the classical 
electromagnetic field), we need to express the atomic nuclei as quantum fields. Then, it is 
necessary to compute the interaction non-perturbatively due to the existence of the bound 
states. Although this is a difficult task to achieve, in this way, we can go beyond the quantum 
mechanical treatment with perturbative QED corrections. We believe such a theoretical 
technique opens up a way to study and predict new phenomena. 

To solve the dynamics of the Rigged QED in a non-perturbative manner, we need different 



In Refs. 
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16| . a formalism to simulate the dynamics of the quantum Dirac field has been developed 
but in the absence of the photon field. 
2 Similar use of the word "rigged" is found in a quantum mechanical context as "rigged Hilbert space" 
which means Hilbert space equipped with distribution theory ■ But it should be noted that the concept 
of the rigged Hilbert space is not used in the rigged QED. 
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techniques from those developed for the standard QED. In this paper, as a first step toward 
this issue, we propose a procedure to expand the Dirac field operator by solutions of the 
Dirac equation for electrons in nuclear potential and derive time evolution equations for the 
electron annihilation and creation operators. Similarly, the Schrodinger field operators are 
expanded by nucleus annihilation and creation operators. Then we derive time evolution 
equations for these annihilation and creation operators and discuss how time evolution of 
the operators for physical quantities can be calculated. In the end, we propose a method 
to approximate the operator equations by c-number equations and show some numerical 
results. 

Before ending this section, we note on our notations and conventions which will be used in 



this paper. They mostly follow those of Refs. 



17H19|. We use the Gaussian system of units. 



c denotes the speed of light in vacuum, h the reduced Planck constant and e the electron 
charge magnitude (so that e is positive). We put a hat to indicate a quantum operator 
to distinguish it from a c-number. A dagg superscript is used to express Hermite 

conjugate. A commutator is denoted by square brackets as [A, B] = AB — BA and an anti- 
commutator by curly brackets as {A, B} = AB + BA. We denote the spacetime coordinate 
as x = (x^ 1 ) = (x°, x l ) = (ct, r) where the Greek letter runs from to 3 and the Latin letter 
from 1 to 3. We adopt the convention that repeated Greek indices implies a summation over 
to 3. Other summations are explicitly written. The transformation between contravariant 
and covariant vectors are done by the metric tensor g^ v = diag(l, —1,-1,-1) = g^ u . The 
gamma matrices are denoted by 7 M . 

This paper is organized as follows. In the next section, the framework of the Rigged 
QED which is relevant to the present study is reviewed briefly. In Sec. 1111} we introduce 
annihilation and creation operators for each field operator and describe how we expand the 
field operators. In particular, Sec. IIII Cl discusses how we treat the photon field operator in a 
non-perturbative manner following Ref. 18[. In Sec. IIVI we derive time evolution equations 
for these annihilation and creation operators and discuss how time evolution of the operators 
for physical quantities can be calculated. In Sec. [V] we propose a method to approximate the 
evolution equations of the operators by the evolution equations for the density matrices of 
electrons and atomic nuclei. Under this approximation, we carry out numerical simulation 
of the time evolution of electron charge density of a hydrogen atom. The last section is 
devoted to our conclusion. 
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II. RIGGED QED 



In this section, we briefly review the general setting of Rigged QED 17H19j by showing 
its Lagrangian and equations of motion for the field operators. 

A. Lagrangian 

First, we describe Rigged QED in terms of the Lagrangian. A part of the Rigged QED 
Lagrangian is the ordinary QED Lagrangian. The Lagrangian density operator of QED can 
be written as 

Lqed(x) = -JL-F^(x)F^(x) + L e D e $} ; x) , (1) 

where F fll/ (x) is the electromagnetic field strength tensor which can be expressed by the 
photon field A J 

F^(x) = dpA v (x) - dyA^x). (2) 

We note here that we adopt the Coulomb gauge V • A(x) = in this work. L e is the 
Lagrangian density operator for the electron 



L 



4>, D efl i>j ; x) = c^(x) (ihj^De^x) - m e cj $(x), (3) 



where m e is the electron mass and the Dirac field operator ^(x) represents the electron (and 
positron). The operator with a bar on top is defined by ip(x) = ip' (x)^ . We denote the 
gauge covariant derivative for the electron as 

D ei *(x) = d lt + i-^A li (x), Z e = -l. (4) 

The Rigged QED Lagrangian is this QED Lagrangian "rigged" with the Lagrangian of 
the atomic nuclei which are represented by Schrodinger fields. We denote the Schrodinger 
field operator for the nucleus a by x a (x). This satisfies commutation relations if a is boson 
and anticommutation relations if a is fermion (which in turn is determined by the nuclear 
spin of a). The Lagrangian density operator for the atomic nucleus a can be written as 

La ( jxa, DaOXa, ^«Xa} I x ) = xl( x ) (ihcD a0 (x) + ^-^a( X )^ Xa(%), (5) 
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where m a is the mass of the nucleus a and the gauge covariant derivative of a is 

Da^x)=d, + i-^-A^x) ) (6) 

where Z a is the a's atomic number. Thus, when we have N n types of atomic nuclei in the 

system, 

^RiggedQED^) = --^-F^ u (x)F^(x) + L e ^|^,/) e ^| ] X^j 

N n 

+ 2Z^ a ({Xa, DaOXa, DlXa} ! , (7) 
a=l 

is the Rigged QED Lagrangian density operator. We note that this Lagrangian has U(l) 
gauge symmetry but the Lorentz symmetry is broken by L a . This, however, will not be a 
problem since we are not going to solve the dynamics in a Lorentz covariant way. 



B. Equation of motion 

Here, we show the equations of motion for field operators, ip(x), Xa{x) and A^(x), in 
Rigged QED. They are given by the principle of least action from the Lagrangian density 
operators introduced above. We also define the charge density operator p(x) and the charge 
current density operator j(x). 

We begin with the Dirac field operator ip(x). Since the "rigged" part of the Lagrangian 
L a does not explicitly depend on tp(x), the equation of motion of ip(x) is same as one in the 
ordinary QED. That is 

ih^D e ^(x)ij(x) = m e cip(x), (8) 

if written covariantly, and it can be also expressed in a form 

ih <H^ = ^z^a^x) + a ■ (-ihcV - (Z e e)A(xf) + m e c 2 (^ $(x), (9) 

which has the same form as a Hamiltonian form of the Dirac equation. In our notation of 
the gamma matrices, (3 = 7 and a = 7°7. 

As for the Schrodinger field operator Xa(x), we obtain the equation of motion in the same 
form as the Schrodinger equation of the non-relativistic quantum mechanics as 

d h 2 - 

ih—Xa(x) = -2m D l( x )x°-( x ) + Z a eAQ{x)xa{x). (10) 
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This can also be written as 

= -£ - 2i f fa ■ * - iM) kx) ■ kx) > Ux) 

+Z a eA (x)xa(x), (11) 

— * — * 

where we use the Coulomb gauge condition V ■ A(x) = 0. 

Finally, the equations of motion for the photon field A^(x) have the same form as the 
inhomogeneous Maxwell equations of the classical electrodynamics. In the Coulomb gauge, 
we have 

-V 2 A (x) = 4irp(x), (12) 
- S V.4„W + ( ?S j-V^ W = - jW , (13) 

where p(x) is the charge density operator, j(x) is the charge current density operator. Note 
that in the case of Rigged QED, p(x) and j(x) include the contribution from atomic nuclei 
in addition to that of electrons. These "rigged" charge and current are described below. 

The charge density operator p(x) is the sum of electron charge density operator p e (x) 
and atomic nuclear charge density operator p a {x): 

N n 

p{x) = Pe (x) + J2pa(x), (14) 
a=l 

where 

p e (x) = Z e e4>(x)<y°ti)(x), (15) 
p a {x) = Z a ex{{x)xa{x). (16) 

Similarly, the charge current density operator j(x) is the sum of electron charge current 
density operator j e (x) and atomic nuclear charge current density operator j a {x): 

l(x)=Ux) + J2Ux), (17) 
a=l 

where 

j e (x) = Z e eci>(x)^(x), (18) 

Ja(x) = ^ (^ihxl(x)D a (x)Xa(x) -ih^D a (x)Xa(x)^ -Xaix^j. (19) 



Here, in passing, some notes are in order. First, the equations of continuity hold for each 
species, namely, 

d s. 

— p a (x) + divj a {x) = 0, (20) 

where a = e or a. Second, p a (x) (a = e or a) is connected to position probability density 
operator N a (x) as p a (x) = Z a eN a {x) and j a (x) is connected to velocity density operator 
v a (x) as j a {x) = Z a ev a (x). 

In summary, the rigged charge and the rigged current are 

p(x) = Z e ei>(x)-/°ip(x) + ^2z a exl(x)xa(x), (21) 

a 

j(x) = Z e ecip(x)jip(x) 

+ Uhx\{x)D a {x)xa{x) -ih(p a {x)xa{x)) • Xa{x) J • (22) 

Now, we have a closed set of time evolution equations for field operators i>(x), Xa( x ) 
and A^x): Eqs. (fTTj), flI3]), (ET]> and To solve them, we rewrite them 

using annihilation and creation operators as is done in ordinary QED (or other quantum 
field theories). In the next section, we describe how we expand the field operators by the 
annihilation and creation operators and derive time evolution equations for them. 

III. EXPANSION OF FIELD OPERATORS BY ANNIHILATION AND CRE- 
ATION OPERATORS 

In this section, we introduce annihilation and creation operators for each field operator 
and describe our expansion method. As is explained in the following, definitions for an- 
nihilation and creation operators are different from those used in the conventional QED 
treatment. This reflects our aim to treat dynamics of electrons, photons and atomic nuclei 
in atomic and molecular systems, which are bound states, in a non-perturbative manner. 

A. Electron field 

In QED, the electron is expressed by the Dirac field operator i>{x) and it is usually 
expanded by plane waves, which are solutions for the free Dirac equation. As is well- 
established, this works extremely fine for the mostly considered QED processes which are 
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represented by scattering cross sections. This is because those processes are described by 
a perturbation series to non-interacting system with the interaction mediated by photons 
treated as the perturbation. 

In our case, however, since molecular systems are highly non-perturbative and electrons 
are bounded, the same method is not likely to work. In addition, what we want to know is 
not the cross section, which just describes the transition between infinite past ("in-state") 
and infinite future ("out-state"), but time step-by-step evolution of the systems. To this 
end, we propose an alternative way to expand ^>(x) as 

4>(ct, r) = ^ [e n {t)^\f) + flmi\f)] , (23) 

n=l 

^(ct,f) = ^ [elm% + \r) + £(f)V4 H (*0] , (24) 

n=l 

where e n is the electron annihilation/creation operator, f n (t)/ flit) is the positron 

annihilation/creation operator and ipn{r) (ip„ (r)) are the four-component wave functions 
of the electron (the positron), which are solutions of the time-independent Dirac equation 
for a particle in a nucleus field. Concretely, ip^ (r) can be considered as the ra-th molecular 
orbitals which are obtained by solving the four- component Dirac Hartree-Fock equation for 
the system and Njj is the number of the basis set. 

We note that this expansion of the field operator is similar to the Furry representation 



(picture) of ordinary QED [2ll.l22| in a sense that it uses the solutions of the Dirac equation in 
an external field. However, the annihilation/creation operators in the Furry representation 
do not depend on time (so it can be considered as a variant of the interaction picture) 
whereas those in our expansion carry all the time- dependence of the field operator. In other 
words, we adopt to work in Heisenberg picture. Since we wish to solve the dynamics of 
the system in a non-perturbative manner, the use of interaction picture does not make our 
problem easier. 

Before proceeding further, to make equations below less cluttered, we introduce our 
notation for the annihilation/creation operators and the electron/positron wavef unctions. 
We define 

e n + = e n , (25) 
e n - = fl (26) 
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and 

xj; n+ (r) = ^\t), (27) 
= ^ _) (r). (28) 

Note that the positron creation operator does not have a dagger as superscript in our nota- 
tion. Then the field operator expansion eq. (1231) can be written as 

N D 

4>{ct,r) =J2^^{t)i) n «{r), (29) 

n=l a=± 

and orthonormality of the wavefunctions can be expressed as 

d 3 ripia(r)ip m b(r) = S nm 5 ab . (30) 
The anti-commutation relation can be written by 

&n a i& m b \ = $nm$ab-> (31) 



and anti-commutators of other combinations are zero. Also, for later use, we define the 
electron excitation operator 

S n a m b = $ n a,e m b. (32) 

B. Nucleus field 

We expand the Schrodinger field operator Xa( x ) f° r the atomic nucleus a by the nucleus 
annihilation operator c a i and creation operator c ai as 

N S 

X a (ct, r ) = d ai(t)Xai(r), (33) 
1=1 

N s 

xI(ct,f) = ^cL(t)x: j (f), (34) 
i=i 

where Xai{r) is a set of orthonormal functions 

d^TximXaA^) = hi- (35) 



The commutation relation for c a i and c ai is 



5ij and the anti-commutation relation 



is "fcaijC^I = 5ij. For the other combinations, (anti-) commutators are zero. Which of 
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them is imposed depends on the a's nuclear spin. If a has (half-) integer spin, the (anti- 
commutation relation is imposed. 

As is done for the electron case, we define the nucleus excitation operator 

C a ij = ^ ai C a j, (36) 

for the later convenience. 



C. Photon field 

In QED, the photon is expressed by a vector field operator A^x) with U(l) gauge symme- 
try. In the standard QED treatment, the quantization is performed to free electromagnetic 
field and the interaction with the electron is taken into account as a perturbation. Although 
this is not the way we try to solve the dynamics of the system as mentioned previously, 
since we are going to use the expression of the quantized free field as a part of A^(x), we 
begin by showing its expression. We call this A r3id ^(x) to distinguish from the total A^(x). 
In the Coulomb gauge, as is commonly done when quantizing the electromagnetic field (e.g. 
Ref. QflQ), we have 



^rad( C ^5 r) 



0. 



d 3 p 



+tf(p,a)e* k (p, cr)e 



a(p,a)e k (p,a)e 



—icp°t/h^ip-r/h 



icp°t / h g—ip-r/ h 



(37) 



(38) 



where a(p, a) is the annihilation operator of the photon with momentum p and helicity a 
and e is the polarization vector. The photon annihilation/creation operators satisfy the 
commutation relation 



[a(p, a), a\f, a')] = 5{p- p')8 a 



(39) 



and commutators of the other combinations are zero. As for the polarization vector, when 
we take the photon momentum in the polar coordinate as 



^ p° sin 6 cos (f^ 



P 



p° sin 9 sin < 

p° cos 9 
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V 



(40) 



(0 < 6 < 7T, < (f) < 2tt, \p\ = p°), we have 



e(p, a) 



1 

7! 



^cos (f) cos 6 =p z sin 



sin cos 9 ± i cos < 
— sin# 



(41) 



where the double sign corresponds to a = ±1. We note that A° ad (ct, r) and A rad (ct,r^) by 
construction satisfy Maxwell equations f lT2|) and f fl3|) with the right-hand-sides being zero 
{i.e. equations for free fields). 

To solve the dynamics of the system non-perturbatively, we take advantage of the fact 
that the formal solutions of the inhomogeneous Maxwell equations ( I12p and (I13p are known 
fr0m the classical electrodynamics B We shah e m plo y the strategy to express A^(x) by 



such solutions 



18| 



As for the first inhomogeneous Maxwell equation ( fl2l) . since it only contains the scalar 
potential A°(ct,r), its solution is readily written by the sum of the solution of Eq. ( 112"]) 
with the right-hand-side being zero and the particular solution of Eq. ( II 2p . The former is 
A|? ad (d;, f) which is zero (Eq. ( l3~Tj) ) and the latter is known as the solution of the Poisson 
equation. Thus, 



A (ct,r) 



3 ^p(ct,s) 



\r — s\ 



(42) 



As for the second inhomogeneous Maxwell equation (TT3]) . we can make it simpler by 
decomposing the current j(x) into the transversal component Jt{%) and the longitudinal 
component j'l(x) as 



3{x) 



h(x)+j L [x) 



(43) 



where V ■ = and V x j'l(x) = 0. This transforms Eq. (fl3l) into two equations each 

involving only A°(ct, r) or A(ct, f) 



1 d 2 A V , 4tt^ , , 

d - - i. 

— VA (x) = 4ttj l (x). 



(44) 
(45) 



Since we have already found A (x) as Eq. (142]) . Eq. (14"5]) tells us Jl(x), which in turn 
gives the right-hand-side of Eq. (T44|) by j(r) — jhif)- Then we can obtain the solution of 
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Eq. (|44|) by the sum of the solution of Eq. ( I44j) with the right-hand-side being zero and the 
particular solution of Eq. OH]). The former is A ra d{ct,r) (Eq. (1381) ) which is written by the 
annihilation/creation operators and the latter is known from the classical electromagnetism, 
which we shall call A^ct^r). Written explicitly, 

i A ^= l - (46) 
c J \r — s\ 

where 

u = t- . (47) 

In summary, the vector potential A(ct,r) is 

l(ct, r) = l rad (ct, f) + A A (ct, r), (48) 
which are given by Eqs. fl38|) . (fISj) and (j4"Tj) . 



D. Charge density and charge current density 

In this subsection we express the charge density operator p(x) (Eq. (12T]) ) and the charge 
current density operator j(x) (Eq. fT22|) ) in terms of annihilation and creation operators (or 
excitation operators) which are defined in previous subsections. 

The expression for the electronic charge density operator can be derived by substituting 
Eq. ( 1291) into Eq. (I15j) . and the atomic nucleus charge density operator by Eq. ( )33l) into 
Eq. ( fT6l) . They give 

Pe(x) = PpV( f )4V' ( 49 ) 

p,q=l c,d=± 
V S 

Pa(^) = ^ Paij{r)C aij , (50) 

where 

p pC(?d (r) = (Z e e)^ c (f)^(r), (51) 
Paii(r) = (Z a e)x* ai (f)Xaj(f)- (52) 
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The expression for the charge current density operator can be derived in a similar manner. 
From Eq. (fT8j) . the electronic charge current density operator is 

N D 

= E E i^(f)4v. ( 53 ) 

p,q=l c,d=± 

where 

j k pCqd (r) = Z e ec [4< (r) 7 VVy (**)] • (54) 
From Eq. (fT9|) . the atomic nucleus charge current density operator is 

N S 



E { Jaiji^Caij — — —Paijif) (^fad^aij + C^A^C^ 1 , (55) 



ij=l 

where 

4-(f) = -(^e)^- ta(f)V*x*(f) - (V^r)) Xa ,(r)} . (56) 

Note that in our notation the spatial components of the covariant derivative are D(x) = 
— V + iS-A(x). Also, we have used the fact that c a j commutes with A k &A but not with A k A . 

Therefore, the total charge density operator p(x) and the charge current density operator 
j(x) are 

N D N n N S 

p( x ) = E E p^a^p^ + E E p*n(?)c«v> ( 57 ) 

p,q=l c,d=± a=l i,j=l 

N D 

"j k ( x ) = E E ^v( r ~*)4v 

p,q=l c,d=± 

N n N S , (7 \ / \\ 

+ E E -4 - -^v p «^ ( Ak ^ + ■ ^ 

a=l i,j=l v y 

Now, we can rewrite the scalar potential A°(ct,f) using the excitation operators. Substi- 
tuting Eq. (J37D into Eq. (J42]) gives 



Af D iV n N s 

A (ct,r) = E E VpcqiiflEpeqd + E E V aij(^aij, (59) 
p,q=l c,d=± a=l i,j=l 



where we define integrals 



V pCqd (R) = I d 3 s^^-, (60) 
Is — R\ 



V aij (R) ee I d's^M-. (61) 
s — R\ 
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Then, we can express ji(x) using the excitation operators via Eq. ( )45l) . Substituting Eq. ( )59l) 
into Eq. ( H5l) gives 

iiw = - E E ^vM^ - E E ^(o^. (62) 

p,q=l c,d=± a=l i,j=l 

where we define integrals 

Finally, we can obtain jr{x) from Eqs. (158]) and fl62|) . and in turn AA(ct, f) using Eq. fj46|) . 
However, the existence of retardation (Eq. (I47p ) in the right-hand-side of Eq. ( H6|) prevent 
us from making the expression of AA(ct,f) simpler. In other words, ^(ct, r) is determined 
by 3t at earlier times than t, whose expression contains at these times. Thus, ApXct^r) 
is computed step-by-step in time using all the information of itself at earlier times. 



IV. TIME EVOLUTION OF ANNIHILATION AND CREATION OPERATORS 

In this section, we derive time evolution equations for annihilation operators of the elec- 
tron e p c and the nucleus c a j. Those for the creation operators can be obtained by taking 
Hermite conjugate. We note that in our formalism, the photon annihilation operator a(p, a) 
does not depend on time as is in the standard QED treatment. We also discuss the time 
evolution of the excitation operators and physical quantities. 



A. Electron annihilation operator 

We first derive the time derivative of the electron annihilation operator e p c This is 
obtained by substituting the expansion f l2U|) into multiplying by ijr p c{f) and integrating 
over r. Then, the orthonormality of the wavefunctions (Eq. (15Uj) ) yields 

iK—r^j- = {jlp c q d + ^2p c q d + ^3p c q d + ^Ap c q d ^ &q d i (65) 

<j=l d=± 
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where 



lp c g d = I d 3 fiplc(f)a- i—ihcV\ip g d(r), 



l 2p c q d 



l 3p c q d 



d 3 f^Uf)a • ( -{Z e e)A{x) ) Vy(r), 



triple (r) (m e c 2 )(3tp q d (f) , 
d 3 f {r) (Z e e) A Q (x) ip q d (f) . 



(66) 
(67) 
(68) 
(69) 



Ii p c q d and I^pcqd are contributions from kinetic energy and mass energy respectively. We 
define the electron kinetic energy integral 



1 pcqd 



—ihc I d 3 r '0pc(r)7 O 7 • V^y(r), 



and the electron mass energy integral 

M pCq d = m e c 2 I d 3 f-ippc(r)'j -ip q d( y f). 



Then, we have i\ p c q d = T pCq d and I 3p c q d = M pCq d. 
i^pc q d can be rewritten using Eq. (1591 as 



l^pcqd 



d 3 rppc q d(f)A (x) 

N D N n N S 

E E wv e s f )£r° 8 f + E E (pviw.)£ 

r,s=l e,/=± 0=1 i,j=l 



aij, 



where we define two types of four-center integral as 



(pY|rV) = (Z e ej 2 I d 3 rd 3 s^ pC (r)ilj q d(r)— — ^(^(s) 



/,t 



r — s\ 



ip c q d \iaja) = (Z e e)(Z a e) \ d 3 r d 3 s^\,{r)^ q d{f) ^^ X*ai(s)Xaj(s)- 



We can write -^pV using Eq. ( 1511) as 



^2p c q d ' 11 



— d rj p c q d(r) ■ (A rad (x) + A A (x) 



Here, the term including A ra d(x) can be written, by substituting Eq. ( 1581) . as 
1 



d rj p c q d(f) ■ A iad (x) 



1 VAnWc 



x 



F pcgd (p) • e(p, a)e- icp0 */ ri a(p, a) + • e^p, a)e icpH ' n a\p, a) 



(70) 



(71) 



(72) 
(73) 

(74) 
(75) 

(76) 



(77) 
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where we define the integral which is the Fourier transformation of the function for the 
electron charge current density (Eq. fl54|) ) 



F k rq M = J d 3 rj k pcqd (r)e^ h . (78) 
The term including Aa(x) can be written, by substituting Eq. ( 146|) . as 

- \ j d\% cqd (rl ■ Mx) = -I J ^^ ^^^'^ (79) 



Further simplification is not possible due to the retardation u = t — in j T . 
Putting these all together, we obtain 

de Nd Nd 
ih ~ftj- = z2z2( T P c q d + M p c qd )e q d + Y (P c q d \ resf )£r*sfe q d 

q=l d=± q,r,s=l d,e,f=± 

N D N„ N S 

+ Y Y Y Y (p C 1 d \ i aja)Cai j e q d 
q=l d=± a=l i,j=l 

- ^2^2^ / rfrrfs yrrgi V 

5=1 <f=± J I 1 I 

1 -\/47r^ 2 c /" c? 3 p 



Fpc,d(p) • e(p, a)e- icp ° t/;i a(p, a)e ? d + F p c qd (-p) ■ e*(p, a)e icp ° t/h a\p, a)eJ . (80) 
The Born-Oppenheimer approximation version of this equation is derived in the Appendix 1X1 

B. Nucleus annihilation operator 

We next derive the time derivative of the nucleus annihilation operator c a j. Similarly to 
the electron case, this is obtained by substituting the expansion into ( TIT]) , multiplying 
by xlii^) an d integrating over r. Then, the orthonormality of the expansion functions 
(Eq. d35])) yields 

N S 

dt' = 



Ilaij + ^2aij + ^aij + liaij^J C a j, (81) 
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where 



*-2 r 

Lij = -^-J d'fxU^Xa^, (82) 
4* = ^ / ■ (x:,(rlV Xa ,(r)) , (83) 

Tfl a C J \ / 

4* = / d 8 *^*) •l(x)x: j (r)Xa J (r), (84) 



2m a c 2 

h mj = (Z a e) J d 3 fA (x) X * ai (f f )Xa J (^. (85) 

/icy is the contribution from nucleus kinetic energy. Defining the nucleus kinetic energy 
integral by 

Tan = / rf 3 r X :,(r)V 2 Xa,(r), (86) 

we have I laij = T aij . 

I iai j can be rewritten using Eq. fl59l) . 



^ = J ^(^w < 87 > 
= E(pViw«)^+^5](u«iyi)4i, (88) 

p,q=l c,d=± 6=1 fc,Z=l 

where we use Eq. ([75]) and define another four-center integral 

(U7a|W = (Z a e)(Z b e) J d'rd'sxU^XaA^r^rXtki^Xbli^- (89) 
iiaij can be written as 

Uaij = ~\j d ^3aij{r) ■ (A^ d (x) + A A (x)} , (90) 

where we integrate by parts, use the Coulomb gauge condition and Eq. ( )56|) . We note that 
this has the same form as Eq. fl76|) but j p c q d (r) replaced by j ai j (r) . Then the term including 
A-ad(^) should have the right-hand-side of Eq. ( J771) but F pCq d{p) (Eq. ( J78l) ) replaced by 
F a y(p) where 

*S,(p) = I d^rj k aij {r)^\ (91) 



is the Fourier transformation of the function for the nucleus charge current density (Eq. (I56|) ). 
Also, the term including Aa(x) is Eq. (1791 but j p c q d(r) replaced by j a ij{^)- 
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haij can be decomposed into four terms because A ra ^{x) and Aa(x) do not commute. 
Namely, 

t Z n e 



' 3aij 



2m a c 2 



J d 3 f (l rad • I rad + l rad -A A + A A - l rad + A A ■ A A ) Paij (f) • (92) 



Since the last three terms involve Aa{x) which includes the retarded potential, we cannot 
make their expression simpler. The first term can be rewritten using the Fourier transfor- 
mation of the function for the nucleus charge density (Eq. (|52|) ). 

F aij (p) = J d 3 r Pm3 (r)e^\ (93) 

and subsequently defining 

Gaij (p, 0", q, r) = e(p, a) ■ e(q, r) F aij (p + q), (94) 



as 



z a e fx x \ Zat sr f d3 P d3 $ 

J d r [ A - d • A - d ) P ^ T) = A^Hc ? J X 



2m a c 2 J \ ) 4ir 2 m a hc J ^/2f 

G aij {p, a, q, r)e-^ 0+ ^ h a(p, a)a(q, r) + G aiJ (p, a, -q, r)e- ic ^">H{p, a)a\q, r) + 

G aij (-p, a, q, T^-^a^p, a)a{q, r) + G aij (-p, a, -q, T^+^a^p, a)tf(q, r) 
where we used e M (— p, a) = e* fJ, (p, a). 

C. Excitation operators and physical quantities 



,(95) 



As is shown in Eqs. ( 149~|) and ( IBTfl) . operators for physical quantities are expressed by the 
excitation operators S pCq d (Eq. ( 1521 ) and C ai j (Eq. (1561)). (For the field operator expression 
of other physical quantities of our interests, we refer Refs. jl3,[l9|.) Thus, we here show the 



time derivative of £ p c q d and C ai j. 



As for £ p c q d, Eqs. ( 152]) and (165]) (and its Hermite conjugate) give 



ih — J^- = (—eleI r e p oe q d + epcJ g d r ee r e^ , (96) 



r=l e=± 



where we define J p c g d ee I lpCq d + I 2p c q d + I 3pCq d + I 4pCq d and use I T C d = / g d pC . Similarly, Eqs. (j36 
and (1HT1) (and its Hermite conjugate) give 



ih 



dC, 



aij 



dt 

k=l 



\^ak^akiC a j + ^ a JajkC a kj (97) 
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where we define I aij = I laij + J 2o y + I 3o y + I ia ij and use 1]^ = I aji . 

Finally, we describe how we can calculate physical quantities from these equations taking 
the electronic charge density for example. The electronic charge density operator is shown in 
Eq. f )49|) . To make the operator observable, we have to take its expectation value. Since we 
work in the Heisenberg picture as we have mentioned in Sec. IHIl the operator is sandwiched 
by time-independent initial bra and ket. Let |$) be such a ket vector which is constructed 
by multiplying the vacuum |0) by the creation operators at the initial time in an appropriate 
manner for a desired initial condition. Then the electronic charge density p e (ct,r) is 

N D 

p e (ct,r) = ($| : p e (ct,r) :|f)=£ J P P v(r)($| : 4vW : l $ >> ( 98 ) 

p,q=l c,d=± 

where :: denotes the operator in between should be normal ordered. Alternatively, we can 
also say the physical quantity should be defined to have zero vacuum expectation value: 
($| : £ p c q d{t) : |$) = ($>\£ pCq d(t)\$) — (0\£ p c q d(t)\0) . This expectation value can be computed 
using the (anti-)commutation relations after £ pCq d(t) is expressed by the operators at the 
initial time via the evolution equation Eq. f )96|) . To express £ p c q d(t) in terms of the annihi- 
lation/creation operators at the initial time, we may discretize the time variable and solve 
Eqs. f )96|) and fl97|) step by step in time direction. Although this procedure and subsequent 
computation of the expectation value can in principle be carried out, since we are dealing 
with the differential equations of operators, which are not commutative, the algebraic ma- 
nipulation required is hugely demanding. The difficulty would increase very rapidly as the 
time steps grow. We consider such brute force computation is important to obtain as exact 
results as possible and it may not be impossible to do so employing recent developments in 



the field of computer algebra 
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25|. At this stage, however, it is more important to look 



for a method to truck the time evolution of the physical quantities based on Eqs. ( 1961) and 
( 1§T|) with suitable approximations. We will discuss such a method in next section. 



V. APPROXIMATION BY DENSITY MATRIX EQUATIONS 

At the end of the previous section, we have pointed out that the excitation operator 
equations of the Rigged QED, Eqs. fl96|) and fl97|) . are practically impossible to solve by the 
present computational technology. We therefore propose an approximation method, which 
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is described in this section. We also show a result of numerical calculation based on that 
approximation for a hydrogen atom. 



A. Density matrix equations 

We begin by introducing density matrices for electrons and atomic nuclei. They are 
defined by the expectation values (Sec. IIV CP of the corresponding excitation operators and 
we denote them as S pCq d and C ai j respectively. Namely, the electron density matrix E pCq d is 

£ pCqd (t) = ($| :£ p c qd (t) : |$), (99) 

and the atomic nucleus density matrix C a ij is 

C aij {t) = ($| : C aij {t) : |$). (100) 

The first step of our approximation is to replace Aa and jx by their expectation values, 
which are denoted by 

A A = ($1 :A A : |$>, (101) 
J T ee ($| : j T : |$). (102) 

Using the density matrices, they are expressed as 

J*(x) = J k (x)-J*(x), (103) 

N D 
p,q=l c,d=± 

+ E E tm^i - —P«v(r) ((Ai A )C aiJ + A k A C atJ ) , (104) 

o=l i,j=l 1 a > 

N D , p N n N s 

= - E E - E E do5) 

p,q=l c,d=± a=l i,j=l 

where (i r fc ad ) = ($|i r fc ad |$), and 

A A (ct, r ) = - / (i 3 s , i ' ; , u = t — ^ !■. (106) 

c J |r — s\ c 

The second step is to replace A Q by its expectation value. This leads to replacement of 
the excitation operators in I 4p c q d and I 4ai j (Eqs. (173"]) and (1851) ) by corresponding density 
matrices. 
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The final step is taking the expectation values of the time evolution equations of the 
excitation operators, Eqs. fl96|) and (1971) . to obtain those of the density matrices. Then, the 
evolution equation of electron density matrix is derived from Eq. ( 19 6 p as 

dS c d Nd 

tfx — = ^ ^ ^ ^ ^E r epc£ r eqd -\- 1rqd r e Sp c r e ^) 5 (107) 

etc 

r=l e=± 

where 1-pCqd = 1.-^pcqd ~\~ l-^pcqd ~\~ L^pCqd ~\~ L^pCqd and 

Xipc q d = Tpcqd, (108) 

l - I d 3 rj p cqd(r) ■ ((i rad (x)) +A A (x)J , (109) 



J -2p c q' 

C _ 

^3p c q d = Mpc q d, (HO) 

N D N n N s 

X 4pV = Yl SbVlrV^Z + ^^^Vliaja)^. (HI) 

r,s=l e,/=± a=l i,j=l 

Similarly, the evolution equation of atomic nucleus density matrix is derived from Eq. ((9 

as 

N S 

oft 



^ ^ ^ ( -E-akiCakj ~\~ ^ajk^aik) (H^) 

k=l 



where we define l aij = l laij + l 2aij + l 3aij + l 4aij and 

(113) 
(114) 

J d 3 f ((l rad • l rad ) + 2(l rad ) • Xi + Xi • Aa) (H5) 

Z^aij = ^ ^ (pV|*oJo)^gd + 5353(iajo|M6)Cwfci. (116) 

p,q=l c,d=± 6=1 k,l=l 

Eqs. (11071) and f 1 1 1 2 [) . together with Eqs. (I103I) - (I106I) . form a closed set of time evolution 
equations of the density matrices. Since they are c-number quantities, they can be solved 
by a straightforward manner. 



7" 

J -laij 


— T ■ 

- 1 at] i 


-£-2aij 




Z a e 
2m a c 2 



B. Numerical calculation 



We here show the result of numerical computation of solving the equations derived in 
Sec. IV Al We adopt to make the equations simpler by employing further approximations. One 
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— * 

is neglecting Aa contribution from the vector potential to avoid heavy numerical integration 
involving the retarded potential. Another is the Born-Oppenheimer approximation (Sec. [A} 
to eliminate the degree of freedom of the atomic nucleus density matrix. Then, we need to 
only solve the equation for the electron density matrix Eq. f ll07j) with T pCq d being 

X pCq d = h pCgd + (pVI^)^/ / d 3 rj p c gd (r) ■ (A iad (x)), (117) 

r,s=le,/=± *' 

where the first term on the right-hand side is defined by Eq. (1A8|) and the third term can 
be written as Eq. (ITT)) with a and a) replaced by (a) and (a?) respectively. 

Now, we describe our computational setups for solving time evolution of electron charge 
density of a hydrogen atom. The orbital functions which are used to expand the electron 
field operator (Sec. IIII Aj) are obtained by the publicly available DIRAC 10 code j^j. The 



Dirac-Coulomb Hamiltonian and the STO-3G basis set are used. Note that in this basis 
set, the electron density matrix is 4 x 4 matrix whose components denote electron (1 + ), its 
Kramers partner (1 + ), positron (1~) and its Kramers partner (1~). The initial state for 
the electron is taken to be the ground state. In terms of the electron density matrix, at 
the initial time, £1+1+ = 1 and the other components are set to be zero. As for the initial 
state for the photon, we perform calculation with and without initial photon field. When 
we include initial photon, we assume coherent state for the initial state so that (a) and 
(at) have non-zero values. We take the coherent state to be that of a single mode with the 
eigenvalue equals to 1. The mode is chosen to be in the x-direction (8 = 7r/2 and = 0) 
and the polarization to be a = +1 (see Sec. IIII CI) . The momentum p° is taken to be 1, 10 
and 20. We work in the atomic units so that m e = e = h = 1 and c = 137.035999679. The 
la.u. of time corresponds to 2.419 x 10~ 17 s or 24.19 as. 

As is described in Sec. lIV C] the electronic charge density p e (ct, r) is calculated as Eq. (I9"S1) . 
In our approximation by the density matrix, 

N D 

p e (ct,r) = PpV^^pvW- ( 118 ) 

p,q=l c,d=± 

We compute this quantity at (x, y, z) = (1, 0, 0) and plot the time evolution in Figs. HH 
Fig. [T] shows the result with no photon in the initial bra and ket. Figs. [21 E] and |4] show the 
result with the initial photon bra and ket taken as a coherent state (as describe above in 
detail) with p° — 1, 10 and 20 respectively. 
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The common feature of these results is the oscillatory behavior. From the visual inspec- 
tion of the figures, there are two types of oscillations. First, we see very rapid oscillations 
with the period of about 1.7 x 10~ 4 a.u. in all of Figs. HH Second, as is seen in Figs. [3] 
and HI there are oscillations with longer period which modulate the rapid oscillations. They 
have the period of roughly 4.6 x 10~ 3 a.u and 2.3 x 10~ 3 a.u in Figs. |3J and H] respectively. 
These numbers suggest the origins of the oscillations. The period of rapid ones is very 
close to the period which is determined from mass scales of electron and positron, that is, 
27r/(2m e c 2 ) = 1.67 x 10~ 4 . The longer periods seen in Figs. |3]andH]are very close to the 
period which is determined from the initial photon momentum, 2ir/(p c), giving 4.59 x 10 -3 
and 2.29 x 10 -3 for p° = 10 and 20 respectively. The interpretation of the latter oscillations 
is that they are caused by the initial photon state which operates as the external oscillat- 
ing electromagnetic field. Technically, they originate from the third term of Eq. fll 17j) and 
the time scale of the oscillations is easily read off from Eq. ( I77|) . As for the former rapid 
oscillations, since it has the period of twice the mass of electron (or positron), it can be 
interpreted as the fluctuations originated from virtual electron-positron pair creations. The 
numerical origin is of course the mass term, Eq. ( ITT]) , in the first term of Eq. ( 11171) . To 
see this electron-positron oscillations more explicitly, in Fig. [51 we plot the time evolution 
of (1 + , l~)-component of the electron density matrix in the case of no photon in the initial 
state. (Note that 1 + denotes electron and 1~ positron as explained earlier.) We see the 
oscillations with the period same as those in Fig. [TJ Although some of the other components 
exhibit similar oscillatory behavior, their amplitudes are much smaller and almost all of the 
contribution to the oscillations come from the (1 + , l~)-component. 

VI. CONCLUSION 

In this paper, we have discussed how we formulate time evolution of physical quantities in 
the framework of the Rigged QED treating non-perturbatively the interactions among elec- 
trons, atomic nuclei and photons. We have defined the time-dependent annihilation/creation 
operators for the electron and the atomic nucleus from corresponding field operators by ex- 
panding them by appropriate functions of spatial coordinates. The photon fields have been 
shown to expressed by the photon annihilation/creation operators of the free photon field 
and those of the electron and the atomic nucleus. This enabled us to include the interactions 
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mediated by the photon field in a non-perturbative manner. We then have derived the time 
evolution equations for the excitation operators of the electron and the atomic nucleus and 
sketched how physical quantities can be computed from them. 

We have pointed out that the last parts of the procedure are very computationally de- 
manding. Solving the coupled evolution equations for the electron and nucleus excitation 
operators by a finite difference method in time requires sophisticated coding technique due 
to non-commutativity of the operators. The existence of the retarded potential makes the 
time evolution method more difficult. Also, since the coefficients of the equations involve 
spatial integration of up to six dimensions, their computation is very time-consuming. After 
that, similar difficulty of non-commutative algebra lies in computing expectation values of 
the time-evolved excitation operators to obtain time-evolved physical quantities. This is 
somewhat alleviated by the use of the Wick's theorem but it still requires large compu- 
tational resource even with small number of time steps. It is true that such brute force 
evaluation can be carried out for a few time steps for a small system and would be feasible 
in future for longer time scale and larger systems considering the recent development in 
computer science and technology. Also, it is certainly meaningful to obtain accurate results 
in such a way. However, at this beginning stage, we consider it is more important to look 
for approximation methods to truck the time evolution of the physical quantities without 
too much computational time. 

Therefore, we have proposed a method to approximate the time evolution equations of 
the operators by replacing some operators with their expectation value in the evolution 
equations. In other words, in an exact sense, we have to compute expectation value after 
the operators are evolved, but, in our approximation, some operators are replaced by their 
expectation values and expectation values are evolved. Although this method has room for 
improvement, we consider this is a good point to start. In fact, the numerical result of the 
time evolution of charge density of a hydrogen atom exhibits the oscillatory feature which 
is considered to originate from the electron-positron pair creations. This shows that our 
simplified way of computation can reproduce one of the most notable features of QED. Even 
within the framework of this approximation, there are a lot more works to be done. In our 
near future work, we will remove the Born-Oppenheimer approximation and will include the 
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full effect of vector potential. 
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Appendix A: Born-Oppenheimer approximation 

In this section, we show the time evolution equation of the operators for the electron under 
the Born-Oppenheimer (BO) approximation. Namely, we fix positions of the nuclei. Since 
the Rigged QED is proposed to include nucleus motion, this approximation is somewhat 
contradictory. However, there are many phenomena which can be described under the BO 
approximation. Moreover, since it eliminates the nucleus degree of freedom, the equations 
become extremely simpler and more accurate results can be obtained with less difficulty 
compared with the non-BO case. Thus, we consider it would be useful to present here the 
BO approximated version of the evolution equations. 

Under the BO approximation, the atomic nuclear charge density operator p a {x) (Eq. ( 1TB]) ) 
is given by 

p a (x) = Z a e5(f-R a ), (Al) 

where R a is the position for the nucleus a, R a = R ai © R a2 © • • • © R an specified by the 
generic direct sum of c-numbered n a vectors that are clamped in space if any, with which 
mandatory manipulation for a function / of R a should read f(R a ) = Ylk=i f(Ra k ) as made 
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obvious, and the atomic nuclear charge current density operator j a (x) = 0. 

This approximates the operators shown in IIIIDI as follows. Eq. ( |59i) is simplified as 

Z a e 



A (ct,r) =J2J2 V pCqd {r)S p c qd + J2 

p,q=lc,d=± a=l \ r *la\ 



(A2) 



Eq. (J58J) as 



N D 



3 k {x) 



p,q=l c,d=± 



(A3) 



and Eq. (I<32"j) as 



No 
p,<3=l c,d=± 



^ / _,. dSpc q d 



(A4) 



Then A^c^r) can be written as 

1 



p,g=l c,d=± 



]pc q d[s) ij E p c q d(s) d£ p c q d 



|f — s| 



|r — s| (it 



(A5) 



where u = t — 



r— s| 



In addition, Z^c^d (Eq. ( J73l ) is simplified as 

I, 



= Y (P C( l d \ resf )^sf + ^2(Z a e)V p c q d(R a ). 

r,s=l e,/=± o=l 



(A6) 



Therefore, the BO approximation version of Eq. flHUl) can be written as 

N D N D 

^2^2h p c qd e q d + (p c q d \r e s f )£ r e sf e q d 



dt 



q=l d=± 



q,r,s=l d,e,f=± 



1 Nd r 



3p°qd{r) • 3re s f{s) £ , s . j p oqd(r) ■ E r e s f(s) dS, 



q,r,s=l d,e,f=± 



\f — s\ 



-£ r e s f(u) + 



|r — s| 



m > e„d 



d 3 p 



x 



F pCgd (p) • e(p, a)e-^ n a(p, a)e qd + F pCqd (-p) ■ e*(p, a)e^ h a\p, a)e qA 
where we define 



h p c q d T p c q d + M p c q d + ^ ^j yZ a t) VpCqd^Rg) . 



(A7) 



(A8) 



a=l 



Under the BO approximation, this is the only equation which governs the time evolution of 
the system. 
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FIG. 1: The time evolution of charge density of the hydrogen atom at (x,y,z) = (1,0,0). The 
variation from the initial value is plotted. The atomic units are used. There is no photon in the 
initial state. 
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FIG. 2: Similar to Fig. Q]but the initial photon state with p° = 1. 
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FIG. 3: Similar to Fig. [I] but the initial photon state with p° = 10. 
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FIG. 4: Similar to Fig. [T] but the initial photon state with p° = 20. 
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FIG. 5: The time evolution of the (1 + , 1 )-component of the electron density matrix for no photon 
in the initial state. 
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